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ABSTRACT 


An.  introductory  discussion  of  the  mathematics  behind  the  direc¬ 
tional  analysis  of  ocean  waves  is  presented.  There  is  sufficient 
detail  for  a  reader  interested  in  applying  the  methods;  further,  the 
report  can  serve  as  an  entry  into  the  theory.  The  presentation  is 
basically  tutorial  but  does  require  a  reasonably  advanced  mathe¬ 
matical  background.  Results  of  a  program  for  the  measurement  of 
directional  ocean  wave  bottom  pressure  spectra  are  Included  as  an 
appendix.  This  second  edition  makes  corrections  to  the  first  and 
adds  some  details  of  an  iterative  directional  analysis  method. 
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1 .  INTRODUCTION 


The  report  presents  an  introductory  discussion  of  the  mathematics 
pertaining  to  the  directional  analysis  of  ocean  waves.  The  presenta¬ 
tion  is  tutorial  in  form  but  does  require  a  reasonably  complete  mathe¬ 
matical  background;  a  background  equivalent  to  that  required  in  reading 
Kinsman's  textbook  Wind  Waves  (1965). 

The  level  of  the  presentation  is  moderate  at  the  beginning.  The 
level  picks  up  rapidly  toward  the  middle  but  there  should  be  sufficient 
detail  and  redundancy  in  the  mathematics  to  allow  the  reader  to  follow 
the  development  without  having  to  rediscover  too  many  omitted  steps. 

It  is  in  this  sense  that  the  report  is  tutorial.  In  some  places  the 
mathematical  development  is  intuitive  rather  than  rigorous.  This  is 
deliberate  in  order  to  provide  insight  and  understanding.  In  most  such 
cases,  references  to  rigorous  reports  are  given. 

The  development  is  reasonably  detailed  so  that  the  interested 
reader  may  apply  the  methods  presented  and  use  the  report  as  an  entry 
point  into  the  rigorous  theory  of  the  directional  analysis  of  ocean 
waves.  In  this  respect,  if  the  report  serves  as  a  bridge  across  the 
gap  between  a  handbook  and  a  rigorous  and  sparse  theory  on  the  subject 
then  the  objective  of  the  report  will  have  been  fulfilled. 

The  report  first  presents  an  intuitive  development  of  a  sea  sur¬ 
face  model  that  assumes  the  sea  surface  to  be  a  two-dimensional  random 
process  definable  in  terms  of  a  directional  power  spectrum.  A  discus¬ 
sion  of  the  space  and  time  covariance  function  and  its  relationship  to 
the  directional  power  spectrum  follows.  Both  one-  and  two-sided  power 
spectra  are  discussed;  however,  the  main  development  is  in  terms  of 
the  two-sided  spectrum.  Next,  the  relationship  between  the  power  and 
cross  power  spectrum  for  two  fixed  locations  and  the  sea  surface 
directional  spectrum  is  developed.  Explicit  relationships  for  the 
special  cases  of  an  isotropic  sea  and  a  single  wave  of  a  given  direc¬ 
tion  and  frequency  are  then  obtained.  The  related  topic  of  the  direc¬ 
tional  resolving  power  of  an  array  of  wave  transducers  is  then 
presented. 

Using  the  preliminary  developments  as  a  basis,  several  methods  for 
the  directional  analysis  of  ocean  waves  based  on  the  information 
obtainable1 from  an  array  of  wave  transducers  are  presented.  The  methods 
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are  basically  a  direction  finder  technique,  a  least  square  single -wave 
train  fit,  and  a  Fourier-Bessel  expansion  fit.  In  conclusion,  a 
generalized  Fourier  expansion  method  is  suggested.  Extensive  results 
of  the  application  of  the  least  square  single -wave  train  fit  are  pre¬ 
sented  in  Appendix  A.  Appendix  B  is  a  FORTRAN  II  listing  of  a  program 
for  this  analysis. 


2 .  WAVE  MODELS 


In  its  simplest  form  an  ocean  wave  can  be  thought  of  as  a  single 
frequency,  sinusoidal,  infinitely  long  crested  wave  of  length  X,  moving 
in  time  over  the  ocean  surface  from  a  given  direction  0.  Such  a  wave 
is  illustrated  in  Figure  1. 


Spatial  Frequency  T1(u,v,t0)  Spatial  Frequency  Along 


FIGURE  1.  SIMPLE  OCEAN  WAVE 
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Assume  that  the  wave  is  frozen  in  time  over  the  surface  (the  x,y 
spacial  plane).  The  coordinates  (u,v)  are  a  6  degree  rotation  of  the 
(x,y)  coordinates.  The  positive  u  axis  lies  along  the  direction  from 
which  the  wave  is  traveling.  The  wave  surface  n(u,v),  shown  frozen  in 
time  in  Figure  1  can  be  described  mathematically  by 

n(u,v)  =  cos(2ttKu  +  2ir<j>)  (2.1) 

where  K  =  1/X  is  the  wave  number  of  spacial  frequency  in  cycles  per 
unit  length  along  the  u  axis,  and  2tt$  is  a  spacial  phase  shift. 

To  make  the  wave  move  in  time  across  the  spacial  plane  with  a 
time  frequency  f  =  1/p,  where  p  is  the  wave  period,  it  is  necessary  to 
add  a  time  part  to  the  argument  of  the  cosine  function  in  the  model 
above.  The  time  part  is  a  phase  shift  dependent  only  upon  time.  As 
time  passes,  the  time  part  changes  causing  the  cosine  wave  to  move 
across  the  (u,v)  plane,  in  this  case  the  ocean  surface.  Adding  the  time 
part  we  get  (where  2mp  is  a  fixed  time  phase  shift) 

n(u,v,t)  =  cos(2ttKu  +  2n<p  +  2ir  ft  +  2ir\p)  . 

If  we  combine  the  effect  of  the  <)>  and  \p  phase  shifts  as  a  =  <p  +  \p, 
we  get 

n(u,v,t)  =  cos  (2rr (Ku  +  ft  +  a))  (2.2) 

as  a  simple  model  of  a  sinusoidal  wave  moving  in  time  over  the  ocean 
surface . 

Since  the  coordinates  (u,v)  are  a  rotation  of  the  coordinates  (x,y) 
through  an  angle  of  0  degrees,  we  know 

u  =  x  cos  0  +  y  sin  0 
v  =  -x  sin  0  +  y  cos  0. 

Using  the  above  relations,  and  letting  £  =  k  cos  0  and  m  =  K  sin  0 
be  the  spacial  frequencies  along  the  x  and  y  axes,  respectively,  we 
have 

n(x,y,t)  =  A  cos(2tt(£x  +  my  +  ft  +  a))  (2.3) 

as  a  model  for  a  wave  of  height  2A  moving  from  a  direction 
0  =  arctan  (m/£) 

with  a  phase  shift  of  2ira.  A  wave  crest  of  such  a  wave  system  is  infi¬ 
nite  in  length.  A  crest  occurs  at  a  set  of  points  (x,y,t)  which  satisfy 
the  relation 


£x  +  my  +  ft  =a  constant  =  (n  -  a) 
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where  n  =  0,  -1,  +1,  -2,  +2,  . . .  .  Each  value  of  the  index  n  relates  to 
a  particular  crest.  The  intersections  of  the  crests  with  the  x  and  y 
axes  move  along  the  respective  axes  with  time  velocities  Vx  =  -f/£ 
and  Vy  =  -f/m.  This  follows  from  the  differential  expressions 


-tfy  -j  <2-4> 

(2.5) 

t  >n  *v*  J  ^ 


obtained  from  the  wave  crest  relationship  given  above. 

From  Euler’s  equation  we  know  that  cos  y  =  (Exp(iy)  +  Exp(-iy))/2. 
If  we  consider  y  as  2tt(£x  +  my  +  ft  +  a)  we  can  write 

n(x,y,t)  =  1/2A  Exp(i27r(ix  +  my  +  ft  +  a)  + 

1/2A  Exp  (i2rr(-£x  -  my  -  ft  -  a))  (2.6) 

where  — 00  <  £  <  +  <*>,  — <»  <  m  <  +  00  and  — “  <  f  <  +  «>. 

In  the  above  we  have  introduced  the  notion  of  negative  time  frequencies. 
This  makes  it  possible  to  express  an  elementary  wave  in  the  mathemati¬ 
cally  convenient  form 

n(x,y,t)  =  a  Exp  (i27r(£x  +  my  +  ft  +  a))  (2.7) 

where  a  =  1/2A.  In  the  real  world  a  complex  wave  of  this  type  implies 
the  existence  of  another  wave  n*(x,y,t)  which  is  the  complex  conjugate 
of  n(x,y,t)  above.  This  complex  conjugate  is  given  by 

n*(x,y,t)  =  a  Exp  (-i2iT(£x  +  my  +  ft  +  a)) 

=  a  Exp  (i2iT.[  (-£)x  +  (-m)y  +  (-f)t  +  (-a)]).  (2.8) 

The  fact  that  negative  frequencies  are  considered  is  explicit  in  the 
above  relation. 


A  property  of  the  above  model,  which  will  be  used  later  in  connec¬ 
tion  with  the  directional  analysis  of  waves  from  measurements  obtained 
from  an  array  of  detectors,  is  expressed  by  the  equation  for  the  phase 
difference  of  two  measurements  made  at  two  different  points  in  space 
and  time.  Assume  we  know  the  value  of  q(x,y,t)  at  the  three-dimensional 
coordinates  (x0,y0,t0)  and  (xQ+X,  y0+Y,  tQ+T) ,  where  X,  Y,  and  T 
are  constants.  The  phases  at  the  two  points  are  given  by- 
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(2.9) 


<t>(xo»yo»to)  =  £xo  +  myo  +  fto  +  “• 

<p(xQ  +  X,  yo  +  Y,  t0  +  T)  £(xo  +  X)  +  m(yo  +  Y)  +  f(tQ  +  T)  +  a 

(2.10) 


This  gives  a  phase  difference  of 

A<j>  =  (£X  +  mY  +  fT).  (2.11) 

To  obtain  a  more  complicated  wave  system  consisting  of  many  waves 
of  various  frequencies  and  directions,  we  can  linearly  superimpose  (add 
up)  many  waves  of  the  form  given  above.  If  we  do  this,  we  can  write 


(2.12) 

*■1 


For  this  wave  system  to  be  real,  the  terms  must  occur  in  complex  con¬ 
jugate  pairs  as  indicated  above. 

For  completeness,  consider  a  model  for  an  infinite  but  countable 
number  of  distinct  (discrete)  waves  and  write 


=  +  .  (2.i3) 

*!•  I 

Again  the  terms  must  occur  in  complex  conjugate  pairs  for  the  wave 
system  to  be  real.  This  will  be  assumed  to  be  the  case  in  future 
discussions . 

A  model  for  a  wave  system  in  the  case  where  energy  exists  for  con¬ 
tinuous  intervals  of  frequency  and  direction  should  be  considered.  In 
particular,  consider  the  general  case  of  continuous  direction  from 
0  to  2tr  radians  and  continuous  frequency  in  the  interval  (-fn,  +  fn) » 
or  even  the  interval  (-“,  ») .  In  theory  the  above  model  does  not  hold 
for  the  continuous  case.  The  power  spectrum  for  the  infinite  but  count¬ 
able  case  would  be  a  set  of  Dirac  delta  functions  of  amplitude  a2 
standing  on  the  points  (£n,  mn,  fn)  of  a  three-dimensional  frequency 
space.  The  continuous  case  produces  a  power  spectrum,  SQ(£,m,f), 
which  is  everywhere  nonnegative  and  in  general  continuous  over  the 
region  of  three-dimensional  frequency  space  where  power  is  assumed  to 
exist.  A  reasonable  model  for  n(x,y,t)  in  the  continuous  case  must  be 
determined.  Consider  a  single  wave  element 

an  Exp(i2ir(£nx  +  n^y  +  fnt  +  an)).  (2.14) 
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The  energy  or  mean  square  in  this  element  is  a£.  Assume  the  ele¬ 
ment  is  a  part  of  a  continuum  of  elements  for  -«  <  f  <  +°°  and  0^6 
<  2tt.  In  this  case  a^  must  be  an  infinitesimal  energy  associated  with 
the  frequency  differential,  df,  and  space  frequency  differentials,  d£, 
and  dm,  which  are  related  to  the  direction  0  of  the  wave  element  as 
before.  Let  the  power  spectrum,  S(£,m,f),  be  defined  with  units  of 
amplitude  squared  and  divided  by  unit  spacial  frequency,  £,  unit 
spacial  frequency,  m,  and  unit  time  frequency,  f.  The  power  spectrum 
is  then  a  spectral  density  value  at  (£,m,f).  In  this  case  we  must  have 
the  infinitesimal  energy,  a^,  defined  by 

s  ^  ^  ^  (2.15) 

The  real  valued,  nonnegative  function  S(£,m,f)  is  a  power  (energy 
density)  spectrum  of  the  standard  type  in  three-dimensional  frequency 
space  (£,m,f).  Intuitively,  we  can  write  an  infinitesimal  wave  ele¬ 
ment  as 

[E%b  (i  M  ( +  XI  ♦  if  *  ,  (2.1.6, 

where  the  positive  square  root  is  assumed.  To  arrive  at  a  model  of 
n(x,y,t)  for  the  continuous  case,  we  need  only  form  a  triple  "sum "of 
the  infinitesimals  or,  to  be  precise,  the  triple  integral  t|(x,y,t)  = 


t  +oc 


(2.17) 


For  different  sets  of  values  of  the  phase  relation  a(£,m,f),  the 
wave  system,  n(x,y,t),  has  a  different  shape,  even  when  S(£,m,f)  is 
fixed.  In  fact  there  is  a  wide  range  of  possible  shapes  of  n(x,y,t) 
for  a  given  S(£,m,f).  The  above  development  is  more  intuitive  than 
mathematically  rigorous.  It  has  been  shown  by  Pierson  (1955  -  pp.  126- 
129)  that  if  2tra(£,m,f)  is  a  random  function  such  that  for  fixed  (£,m,f) 
phase  values  of  the  form  2ttoi  MOD  2tt  between  0  and  2tt,  are  equally 
probable  and  all  phase  values  are  independent,  then  Equation  (2.17) 
represents  an  ensemble  (collection  of  all  probable)  n(x,y,t)  for  a 
given  S(£,m,f).  The  random  process  represented  by  the  ensemble  is  then 
a  stationary  Gaussian  process  indexed  by  the  three  dimensions  (x,y,t). 
Detail  discussions  of  the  above  can  be  found  in  St.  Dennis  and  Pierson 
(1953  -  pp.  289-386)  and  Kinsman  (1965  -  pp.  368-386).  The  fact  that 
a  particular  sea-way  can  be  considered  as  a  realization  of  a  stationary 
three-dimensional  Gaussian  process  has  been  verified.  Refer  to-  Pierson 
and  Marks  (1952) . 

The  model  in  Equation  (2.17)  as  a  stationary  random  process  will  be 
assumed  in  following  discussion. 
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3.  A  DIRECTIONAL  WAVE  SPECTRUM 


The  power  spectrum  S(£,m,f)  is  the  directional  wave  spectrum  of 
n(x,y,t).  If  ri(x,y,t)  is  to  be  real,  every  infinitesimal  of  the  form 
given  in  the  continuous  case  above  presumes  the  existence  of  its  com¬ 
plex  conjugate.  Let  us  consider  the  one-sided  power  spectral  density, 

S ' (£0,m0,fo) ,  of  a  single  real  wave  where  0  £  fQ  <  “.  For  such  a  real 
wave  element  of  length  XQ,  from  a  direction  0Q,  0  <_  0Q  <  2ir,  the  value 

S'  (£0>m0.f0)>  =  S(£0,m0,f  )  +  S(-£o>-m0»“fo)>  where  Zo  =  Ko  cos  6o» 
m  =  K  sin0o,  and  K0  =  1/X0.  If  the  above  real  wave  came  from  a 

direction  (2it  -  8)  ,  the  power  density  would  be  S  ao»-mo’fo)  S(£0,-mo,f0) 
+  S (— £q , » ~f q) • 

Figure  2  illustrates  a  real  wave  of  length  XD,  from  a  direction  0Q, 
in  two-dimensional  spacial  frequency  (wave  number)  space. 


m 


If  the  wave  number  relation 

K=  2ltfV8/Ta.HK(aitKli)  (3>1) 

holds,  refer  to  Kinsman  (1965  -  p.  157)  and  Munk  et  al  (1963  -  p.  527), 
where  h  =  water  depth,  g  =  acceleration  of  gravity,  and  K  =  wave  number 
=  1/X ,  X  being  the  wave  length;  then  a  relationship  between  f  and  (£,m) 
is  implied  that  requires  a  wave  frequency  f0  to  have  a  unique  wave  num¬ 
ber  kQ.  From  this  we  have  the  general  one-sided  spectral  form  for  waves 
where  f  =  fQ  of 
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S'  (£,m,f0)  =  /  zero  where  l2  +  m2  j*  KQ2 

(a  power  density  >_  0  for  l2  +  =  KQ^. 

S'(£,m,f0)  thus  defines  power  density  at  f  =  f  for  wave  energy  over 
0  £  0  <  2w.  Figure  3  illustrates  this  case  in  wave  number  space. 


FIGURE  3.  DIRECTIONAL  WAVE  SPECTRUM  AT  A  FIXED  FREQUENCY,  fQ 


We  want  to  estimate  the  shape  of  S’(il,m,f0)  above  the  circle  l2  +  m2  = 
K  2  in  a  directional  wave  train  analysis.  Remember,  the  S'(£,m,f  ) 
a§ove  is  restricted  to  fQ  >_  0  and  is,  in  fact,  equal  to  v* 


where 


t  s *v  s 

»  $  (-  A, -*> ,  - 


(3.2) 


if  n(x,y,t)  is  to  be  real. 


Let  us  see  how  S'  (£,m,f)  might  be  found:  we  have  said  that 
n(x,y,t)  can  be  assumed  to  be  a  stationary  Gaussian  process.  One  char¬ 
acteristic  of  such  a  process  is  that  for  fixed  values  (xQ,  yQ,  tQ)  of 


the  procet 
distribution;  i.e.. 


t0) is  random  variable  with  a  Gaussian 


TVek  (*1  (*• .  & » £ •) *  yl*) =  J  ®  ^  ^  ^  <i 
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where  p  is  the  arithmetic  mean  of  n  and  a  is  the  variance.  Intui¬ 
tively  n  is  as  likely  to  be  positive  as  negative,  so  let  us  assume  that 
Prob  (n(xD,  y0,  t0)  <  0)  =  1/2.  Since  n  is  Gaussian  distributed,  and 
is  thus  symmetric  about  its  mean,  we  have  Prob(n(x0,  yQ,  tQ)  <  p)  =  1/2 
or  that  p  ■  0.  For  ,  we  have  (using  expected  value  notation) 

<r‘-  »..<.)] 

where  wo  are  thinking  of  n  as  a  random  variable. 

A  Gaussian  process  is  completely  defined  statistically  if  we  know 
the  form  of  the  mean 

F(n(x,y,t))  and  the  covariances 

t{[’t t  to(M.t>)]fo(*«X,»VY,t*T)- 

where  X,  Y,  and  T  are  space  and  time  separations,  respectively.  Refer 
to  Parzen  (1962  -  pages  88-89).  We  have  assumed  E(n(x,y,t))  =  p  =  0  and 
that  the  process  is  stationary  (only  weakly  stationary  is  necessary). 
Hence,  by  definition  of  weak  stationarity ,  we  have  for  each  (x,y,t), 
and  yet  independent  of  the  particular  x,y,t  values,  the  covariance  form 

R(UT)s  E[l(u*)  <*T)]  .  (3-3) 

All  of  the  properties  of  the  stationary  Gaussian  process  n(x,y,t)  are 
implicit  in  R(X,Y,T),  just  as  a  knowledge  of  p  and  a ^  for  a  single 
Gaussian  random  variable  completely  defines  such  a  random  variable. 

Here  it  is  important  to  understand  that  we  are  discussing  expected 
values  across  all  possible  realizations  at  a  point  (x,y,t);  i.e., 
across  the  ensemble  of  all  possible  sea  wave  shapes  at  (x,y,t)  for  a 
given  S(£,m,f) . 

There  is  a  simple  and  unique  relationship  between  R(x,y,t)  and 
S(£,m,f).  Consider  a  single  real  wave  element  (from'  Equation  (2.16)  and 
(3.2))  as  a  random  process  and  write  n(x,y,t)  =  [EXP(i2Tr(£xfmy+ft+a)) 

+  ^  ^  f  t  +  XT' . 
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Form  the  covariance  function 


R  =  £[7 (*»*»*)  *  7(«  +  X,ti+Y,  t+T)] 

=s  E  [§*|p  (i  zn  (<  (z*  4  X) + >1  (2  y+ Y) + ■f  (u  +T )  ■*■ *• <*) 
+  £*t»(-i2ir(U***X)*>»»(2^+Y)+,f  (et+T^+a^ 
4*  (i-  Etr  ($X  4  >\  Y  t f  T)) 

4- Ex(p  (-izii  (iX+ T))] 


where  E  is  the  expected  value  over  the  ensemble  for  any  fixed  x,y,t. 
Note  that  -t 


is  a  constant  with  respect  to  the  expected  value. 

Consider  the  following  problem.  Let  u  be  a  random  variable  with 
uniform  probability  density  function 

o  $  U  <  air  r4JUa*4 

else  tolneve 


{: 


Define  a  random  variable  Z  =  e^u.  The  expected  value  of  Z  is  defined 
as 


Considering  the  random  variable  nature  of  the  phase  2ira  as  described 
following  Equation  (2.17)  at  the  end  of  Section  2,  and  applying  the  above 
notion  to  the  cross  product  terms  of  R(X,Y,T)  we  obtain 


R(XYT)  =  [£*j»(i*.ir(IX«-wY*fT))  + 

e  « J.  (_  i  zn  ( I  x  ♦  n r  *  f  T))J  5  (W)  UK.  A  f . 
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For  a  real  wave  element  we  have,  where  S(£,m,f)  =  S(-£,-m,-f),  see 
Equation  (3.2),  R(X,Y,T)  -E^  (i  E1l  Y-  -f  T>  S 


+  E*t»(-iiir(5x+’nY+^T))S(-4-*.,-yJinMAf 


which  is  simply  the  sum  of  the  covariance  functions  of  two  complex 
wave  elements  which  are  conjugate  pairs.  It  also  follows  that 
R(X,Y,T)  is  real  valued. 

Reverting  to  the  complex  wave  element  form,  and  noting  that  the 
expected  ensemble  value  of  cross  products  between  different  wave  ele¬ 
ments  is  zero  in  a  manner  similar  to  the  case  of  cross  products  shown 
above,  we  obtain  the  composite  general  relationship 

Klxyj)-  Jjp *t)S(W)«Uju<»  <3« 

We  have  demonstrated,  but  not  rigorously  proven,  that  the  covari¬ 
ance  function  R(X,Y,T)  is  the  three-dimensional  Fourier  transform  of 
the  directional  power  spectrum  S(£,m,f). 


We  cannot  hope  to  be  able  to  estimate  R(X,Y,T)  for  continuous 
values  of  X,  Y,  and  T.  However,  $here  is  a  way  around  this  problem, 
we  can  write  the  above  as  _ 


(i  ivS (Utt  (Jx*mY)S(a*>.<)aU*,JA4 


(3.5) 


which  is  in  the  form  of  a  single  dimension  (variable  f)  Fourier  trans¬ 
form  of  the  term  in  brackets  [  ]’s.  Note  this  term  is  not  a  function 
of  T.  It  depends  only  on  the  value  of  (X,Y) .  Further,  by  Fourier 
transform  pairs  we  can  write  this  expression  as 


[  ]s  =  (R(xy.T)E*K-'l^^)<n-  (3.6) 

-OP 

In  general,  assuming  that  the  terra  in  [  ]fs  is  complex,  we  can  write 

[C<xy,t) -  Lq  (xxOU  = 

}jMXXT)E*K-iiVT)JlT  •  <3- 
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To  find  [C(X,Y,f)  -  iQ(X,Y,f)]  we  need  only  know  R(X,Y,T)  for  continuous 
T  for  the  given  value  of  X,Y) .  Further,  we  have  just  stated  that 


LcOW)-i©0W)]  =  [  Js  = 

(i  ainix*  *Y))  a 


(3.8) 


This  is  in  the  form  of  a  Fourier  transform;  thus,  we  can  write  from 
transform  pairs 

S(W)  =  jj[c  (XSY,  f)-i  0  (XY  Ol  (-i  (#*>  Mf))<bUY 

(3.9) 

As  has  been  stated,  we  cannot  hope  to  have  a  continuous  set  of  values 
of  (X,Y) .  The  solution  is  to  find  R(X,Y,T)  for  continuous  T  and 
selected  values  of  (X,Y),  and  then  employ  the  above  to  estimate  S(£,m,f) . 
This  is  described  in  the  next  section. 


4.  CROSS  SPECTRAL  MATRIX  OF  AN  ARRAY 


Let  us  look  at  n(x,y,t)  at  two  fixed  points  in  space,  say 
(x0,  yQ)  and  (x^,  y^) .  This  would  give  two  stationary  Gaussian  pro¬ 
cesses  indexed  on  time  alone  because  (x0,  yQ)  and  (x]_,  yj)  are  fixed. 
Thus,  we  may  write 

\M)  =  t (*•**•, <0  n  (*»%£) 

If  X  =  (x^  -  xQ)  and  Y  =  (y^  -  yQ)  Then  we  can  say,  since  n(x,y,t)  is 
assumed  weakly  stationary  (see  Equation  (3.3)),  that 

R(X,Y,T)  =  efr.W-U^T)] 


where  the  expected  value  is  over  the  ensemble  for  some  specific  value 
of  t,  where  -  «  <  t  <  “.  Let  us  extend  this  idea  by  a  change  of  notation 
and  let  N(X,Y,t)  be  a  two-dimensional  (vector)  process,  double-indexed 
on  time;  i.e.,  let  N  be  a  vector  function 


N(XXt)  =  N  W 

for  -oe>  <  t  <  co 


(4.1) 
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We  can  then  write  a  generalized  covariance  function  (assumed  to 
be  finite)  as  the  matrix  equation 


R(T)  =  e.[nV)  N(t+T)J 

_  E.  (*!.  (t)  *  T))  t  (*)  1,  (t  +  T)) 

_E. (7,(*)  HM+T))  E.(l,W 

:R(o,o,t)  R(X,V,T)‘ 

~  Rf-X,-Y,Tl  R(0,0,T\ 


(4.2) 


Now  R(-T)  is 


R(-T)  = 


-T) 

l_R(-XrY.-T) 


R{X,Y,-T) 

R(o,o,-t) 


and  R(-T)  transpose  is 

m 


rt(-t) 


(«.3) 


Rf-X.-Y-T') 

R(X,Y,-T)  R^o,  o, -T~> 

We  have  by  Equation  (3.3)  and  stationarity  that 

R(-X,-Y,-T)  =  E.(»i(*,,S„t+T)  •?  (*.,:(.,*)) 
—  £  (^  fa.,  a. ^  fa/>xM  t * —  R  (X,Y ,T) 


It  then  follows  that 


R  (T)  =  R(-T) 


(4.5) 
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From  Equation  (3.7), 

let  p'(x.Y.t)  =  [C(x.,Y,f)-iQ(x,Y.f)] 

and  we  get 

P*(X,Y,f)  —  (  R(X,Y,T)  E.*\>(-a-MT)<lT  . 

-at 

Now  R(-X,-Y,T)  =  R(X,Y,-T)  by  Equation  (4.5).  Thus,  we  have  from 
Equation  (3.4) 

R  l%yrT)  =jjyEi<^(iair(U  +  ^Y+U-T))5(W)(iUM4f 

or,  by  the  same  procedure,  that  Equation  (3.6)  was  obtained  from 
Equation  (3.4),  we  have 

4b 

P*fX,Y,'f\  =f  R(X,Y,-Tk*|>(iair4T)<!lT  <*-6> 

>«0 

Now,  since  R(X,Y,-T)  is  real  and  Fourier  transform  pairs  are  unique, 
we  must  have  from  Equation  (4.6)  that 

J  R  (X,Y,  -T)  E*|>(-t2f  *t)  at  =  p  (x,Y.  f ) 


a  parallel  form  of  Equation  (3.6).  Therefore,  the  Fourier  transform  of 
R(-X,-Y,T)  =  R(X,Y,-T)  is  the  complex  conjugate  of  the  transform  of 
R(X,Y,T) ;  i.e.,  (see  Equation  (3.7)) 

p*(-x,-Y,  t )  =  [  P*(X,Y.  0]"  =  P  (x.Y.  0 

T 

In  general,  we  have  that  the  Fourier  transform  of  R(T)  =  R  (-T)  is 


~P  (o.o.O  Plx,Y.O 

_P(X,Y.f)  P(o.o.-f) 

Note:  P*(0,0,f)  =  P(0,0,f)  is  real. 


(4.7) 
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By  definition  it  follows  that 
C(X,Y,f)  =  Re[P(X,Y,f) ] 
Q(X,Y,f)  =  Im[P(X,Y,f) ] 
where  -“  <  f  <  °°  . 


The  function  C(X,Y,f)  is  called  the  cospectrum  and  Q(X,Y,f)  is  called 
the  quadrature  spectrum.  Both  are  spectral  density  functions. 

Explicitly,  we  can  think  of  (x0,  yQ)  and  (x^,  y^)  as  being  the 
location  of  elements  of  a  probe-array  with  space  separation  (X,Y) . 

The  first  step  to  find  (or  estimate)  S(£,m,f)  (see  Equation  (3.9)) 
is  to  find  P(X,Y,f)  related  to  a  pair  of  array  elements.  This  is  a 
problem  of  estimating  the  cospectrum  and  quadrature  spectrum  of  a  two- 
dimensional  (vector)  stationary  Gaussian  process.  Goodman  (Mar  1967  - 
Chapter  3)  has  an  excellent  treatment  of  this  subject,  which  we  will 
discuss.  Kinsman  (1965  -  Chapters  7-9)  also  discusses  the  subject. 

The  essence  of  the  problem  is  that  if  P(X,Y,f)  is  continuous  and  negli¬ 
gible  for  | f |  >  fn,  then  R(X,Y,T)  can  be  obtained  by  a  time  average  over 
a  particular  realization  N(X,Y,t)  instead  of  having  to  average  (find  the 
expected  value)  over  the  ensemble.  This  says  that  we  can  find  R(X,Y,T) 
by  obtaining  two  time  series  (realizations)  rQ(t)  and  r^(t)  measured 
over  time  at  only  two  points;  e.g.,  (x0,  yD)  and  (x]_,  yi)  where 
(x^  -  xQ)  =  X  and  (y^  -  yQ)  =  Y.  The  relationship  between  (r  (t) ,  r^(t)) 
and  R(X,Y,T)  -<*>  <  T  <  <*>  is  ° 


Rfx.Y.  T^=  Limit 


/  ^ 

7  )  r0(tV  r,  (t+r)  dt 

-V 


(4.8) 


where  rQ(t)  is  a  realization  measured  over  time  at  (xq,  y  )  and  r-^(t) 
is  measured  at  (x]_,  y^) .  We  can  simplify  the  notation  by  an  expression 
for  a  time  average  given  by 

R(X,Y,T)  =  Roi(T)  =  r0(t)r1(t+T). 

i 

It  follows  from  Equation  (3.7)  that : C(X,Y,f)  =  CQ^(f)  = 

J  Rqi  (T)  cos(2rrfT)dt, 

-03 
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and 


00 

Q(X,Y,f)  =  Qoi(f)  =  J  Ro!(T)  Sin(2TifT)dT.  (4.9) 

-00 


We  also  have  (see  Equations  (4.2)  and  (4.7))  Cqq(:£)  =  C-Q(f)  =  Pg^f)  = 
pll(f) > 

Qoo(f)  -  Qii(f)  -  o; 


Pgl(f)  =  CQ1(f)  -  i  Q01(f).  (4.10) 

The  phase  of  Pg^(f)  is  given  by 

Q0i(f) 

Q01  (f )  =  Arctan  ^ .  (4.11) 


This  is  the  expected  phase  lead  of  the  signal  at  (xQ,  yQ)  over  the 
signal  at  (x^,  y^)  for  f  where  -  <»  <  f  <  ®  . 

For  an  array  of  N  detectors  located  at  (x^  y^  ,  (x2,  y2) ,  (x3,  y3) 
...(xN,  yN)  we  can  find  N2  spectra  [P^]  i  =  1,  N,  j  =  1,  N.  This  gives 
a  unique  spectrum  Pi:L  (Pn  =  ...  =  P  ;  and  since  P±1  =  P*  (see 
Equation  (4.6)  and  4.7))  ^ 

we  have 

(4.12) 


unique  cross  spectra  or  a  total  of  [N(N-l)  +  2]/2  unique  spectra.  Thus, 
for  real  n(x 
us  to  define 
a  cross  spectral  matrix 


»y » t)  we  have  P±j(f)  =  P*j,(f)  and  pi-5^f)  ~  allowing 

the  information  about  P\X,Y,f)  obtainable  from  an  array  by 


C,t(0 

P„W 

•  « 
0 

-  • 


'  C  I  N  (f )  \ 


•  c*«w 


WHtRt  0  it  $  <  OO 


(4.13) 
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This  information  can  be  used  together  with  Equation  (3.9)  to  obtain  an 
approximation  to  S(£.,m,f).  Numerical  details  for  finding  the  spectral 
matrix  are  given  in  Bennett,  et  al  (June  1964). 

It  should  be  pointed  out  that  negative  frequencies  are  still  con¬ 
sidered  in  the  relationships  being  discussed.  We  do  not  know  P*(X,Y,f) 
for  continuous  values  of  (X,Y) .  We  do  know  from  the  spectral  matrix 
the  values  of 

<"  1  *  ••  i*  *•.->» 

where  X.,  =  (x.  -  x. ) 

ij  j  i 


We  also  have  from  Equation  (4.7)  that 

=  T>  (X^.Xj ,  0 

From  Equation  (3.9)  we  get 


(4.14) 


4)  cos  [lit  (8  X  -  ™Y)] 

•oft  -oar 

-  iP?X,Y,<)  siN[an(lXf  mY)j|AxAY 

or,  since  S(jJ,m,  f)  is  real 

$(*.■»». *)  =  \  \  {c Cx  y  *)  + 

—  CO  -*®  ' 

-Q(X.Y.f)siN[an(!X+mY)]}  AXctY  (415) 

Let  us  consider  treating  the  points  of  (X,Y)  where  we  know  C(X,Y,f)  and 
Q(X,Y,f)  as  weighted  Dirac  delta  functions;  e.g.,  at  (X^>  ^12)  we  8et 

C(X,Y,4)=  kltC *)  ^  (Y-  Xxj . 
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Reverting  to  the  C^j ,  Q^j  notation  of  Equation  (4.13),  we  have,  where 
Xj,  =  and  Yjj  =  -Y-h,  Cjj  =  C^.*  and  Q.w  =  -Qi j .  The  numerical 

form  of  Equation  (4*15)  then  Becomes 

S(t,m,£)-blt  Cu(t)  +  tj££  bJC  (<)  costal 

a  =  i  j=i+i  C  ^  J  J 

— •  Q.j(4) 'Sm  [z.^  +  .  <4.i6) 

Choice  of  b^j  values  is  arbitrary.  A  reasonable  choice  is  b..  = 
[N(N-l)  +  l]--*-  for  all  (i,j)  (refer  to  following  section).  We  now^have 
a  basis  for  an  approximation  of  S(£,m,f).  Before  exploiting  this  result, 
we  need  a  few  side  results. 


5.  SPECIAL  CROSS  SPECTRAL  MATRICES 


Assume  that  we  have  a  real  sea  wave  .of  frequency  f  >  0  moving 
from  direction  0O  where  ZQ  -  KgCosQo  and  mQ  =  KQSin  0O,  K0  being  the 
wave  number  from  Equation  (3.1).  We  can  write  the  wave  as 

V  ,1)  —  A  fm.'i i  +  4.  t  4-**})  (5.D 

Since  the  root-mean-square  (rms)  value  of  a  cosine  wave  is  A//27  we  have 
for  the  two-sided  (-»  <  f  <  «)  directional  power  spectrum  of  the  wave  in 
Equation  (5.1) 

~  |y(S-io)i  (>n-7n.)cf(4-{a)  + 

+  i  +  j  i(4  +4-)]  (5.2) 

t  2  2"  JQ 

or  in  polar  form  where  K  =  vZ  +m  ;  0  =  arctan(y) 

X. 

S(K.e.rt=  £[^e-e.)  *(*-*.)♦ 

+  i(e-(s.--n))J(4++.)]i(K-K.) 
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From  Equations  (3.8)  and  (5.2)  we  have  for  a  single  wave  of  frequency 
fD  >  0  from  a  direction  0O  that  the  two-sided 

R’V'l  =  -£•  [e*.p  (i  lit  (S. X..+  m.  Y,))  i (f  - 4.)  + 

+•  EKP  (i (-  1  m  Y;i))  J  (4  ♦  f„)J 
OV*  u>HCR£  P^M  =  C..(4)  —  LQ..($) 


Qj(0=  c/-0  =  v  ^0), 

QjW  =  -QljK)  =  - AlslN(fn(i.  V>n  YAi)) 


(5.3) 


describe  the  elements  for  the  spectral  matrix  of  a  single  wave.  Recall 
that,  in  general,  <M-j(f)  =  cji(f)  and  Qij(f)  =  — Q±-j  )  -  Substituting 

into  Equation  (4.16),  we  get  where  b-j^  =  [N(N-l)+ri  1  =  1/M  ^  ^ 


and  f .  >  0  is  assumed  C/ft.  A  f  i  . 

°  o  (J  =  —  \  l  + 

H  t  r  4^  C 

*2  2,  ( [cos  ( *•  V  j ) cos  ^  (*  V  *  Xi  D  + 

+  S 1 N  (  ZTf  ( l  X-  +  S I N  (  2.Tt  ( 4  *»  %))]  o  s 

S(W^  = 


(5.5) 


The  choice  of  b^j  =  1/M  =  [N(N-1)+1] -1  and  observing  that  ^  K  \  -  HiUzl) , 


t 


gives  the  convenient  result  (Note:  For  -f  <  0  we  would  use  (-£Q,-m0) 
in  place  of  (£0,m0)) 

=  S  (-1,-^0,- (o') 


which  agrees  with  the  values  of  Equation  (5.2)  at  the  points  (£Q,m  ,fQ) 
and  (-£0,-m0,-f0) .  Note  that  there  is  not  general  agreement  elsewhere. 
In  fact,  Equation  (4.16)  may  (and  does)  give  negative  values  for  the 
approximation  of  S(£,m,f).  This  is  a  problem  of  probe  array  design  and 
is  directly  related  to  the  directional  resolving  power  of  a  probe  array. 
This  problem  is  discussed  later  in  another  section. 


Consider  the  case  of  a  single  frequency,  f  >  0,  real  sea  with 
equal  wave  energy  from  all  directions;  i.e.,  isotropic  waves.  In  this 
case,  assuming  the  wave  equation  (Equation  (3.1))  holds,  we  have 


OR  UJH&R.C  J  )  A*J0  ®  =  ARCTAM 

S  (  K  ,0 ,  f )  =  X  i  (k1 -  K.1 )  (f  -f.)  +  l(f +  f.)  ] 


as  the  two-sided  directional  power  spectrum  for  such  an  isotropic  wave. 

Since  £  =  Kcos8  and  m  =  Ksin0,  Equation  (3.8)  can  be  expressed  in 
polar  coordinate  form  as 


if  oo 

T?.(f)  =  j  |S(K,0,4)&xp  (i.  2.TT  K  (X^cos©  -^Y^.SIN  0))KdKd0(5 

i'*  -If  ° 


7) 


Letting 
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and 


we  can  write 


-y-  If  co 

?jW  =  J  (s(K,e,4)EKP(iiwK^jC.*(e-^Kakde  (5.8) 

-K  O 


Using  Equation  (5.6)  for  S(K,0,f)  we  get 

lf(fj  =  %  K.  f  E  ( I  Wos  (e  -  <y) 

=  £  K.  J  Coe  (air  k0Du  cos  (e  -  y )  de 


die 


-T1 
_TT 


+  L  X  K,  jsi  h  (air  K,  D^cos  (e  -  <y)  di  0 


(5.9) 


-It 


Now,  departing  from  the  above  development,  consider  the  following 
integral  where 

z  =  eir  k#d4-3  i  =  <t>ij 

\  cos(»e)  cos  (3i  cos  (e-0)))dle 

-■nr 


Let  <}>  "  9  -  and  we  get 


(p)  cos  (z  cos  $)  dl(J> 
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Expanding  cos(n0  +  ni|/)  we  get 


We  have  (ir  -  i|>)  -  (— tt  -iJ>)  =  2tt  so  that  the  integrands  are  over  2ir 
allowing  us  to  write  the  above  as 


COS 


»  qjf 

t 

—  S I K  4>  J 


COS 


V 

SlN 
-ir 


r\  $  co s(z  cos  4>) i ^ 
n  $  cos  (*  cos  0)  A$  . 


From  Ryzhik  and  Gradshteyn  (1965  -  page  402)  we  have  (since  sine  is  odd 
and  cosine  is  even)  the  result 


^cos  (vie)  cos(i  cos(e- 
-=.  COST)  Jilt  cos  J^OOj  . 


(5.10) 


In  a  similar  way  we  get 


\cos(y>e)  sin(i  cos  (e- i}>))4e 

-tt 

=  cos  SINFUL)  JnU)]  3 


(5.11) 
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where  Jn(Z)  is  a  Bessel  function  of  the  first  kind  of  integer  order  n. 
Returning  to  Equation  (5.9),  the  above  results  give,  for  n  =  0  where 
fQ  >  0,  the  result 


^  (±4J  =  A  Zlt  K.  J.  (2Tt  K.Djfj)  . 

Since  the  above  is  real 

=  ZTr  K.A 

=  ZtrK.ftj.fZ1t  K.D.)  - 

“  0  (5.12) 


c^.O 

QiiUt) 


describe  the  elements  for  the  spectral  matrix  of  a  single  frequency 
isotropic  sea. 

The  above  two  special  cases  for  sea  waves  are  the  extremes  of 
directionality  of  real  sea  waves  of  frequency  fQ  >  0.  These  results 
will  have  important  applications  later. 


6.  A  MEASURE  OF  ARRAY  DIRECTIONAL  RESOLVING  POWER 


From  Equation  (3.9)  we  have  the  Fourier  transform 

S  ( 8 . W .  f )  =  X  0  £*P  (-* w  (j  X+fciY))  ix  dY  (6  1} 

—  CO  -OB 


In  practice  we  use  a  probe  array  with  elements  at  (x^,  y^)...,  (x^,  y^) 
to  obtain  P*(X,Y,f)  at  the  separation  points  (0.0),  (X^,  yi2>  > 

(”Xl2 ’  Y12^,,,»  ^Xk-l,k>  Yk-l,k^’  ^“Xk-l,k>  Yk-l,k^  ’ 

We  do  not  then  know  P*(X,Y,f)  but  the  product 

[P*(X,Y,f)  g(X,Y)]  (6.2) 
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where  g(X,Y)  is  a  set  of  Dirac,  lelta  functions  standing  on  the  separa¬ 
tion  points  of  the  probe  array,  and  zero  elsewhere. 

Thus,  we  have  the  estimate 


S  (9,m,f )  =  jjjV?  x  ,Y  ,4\  3»(x  ,Y)]  EXP  (-  i.  En(ix  +  Y))  Ax  1 Y 

-00-00 


Let 


•O  03 

G  (i ,  w*)  =  J[^(x  ,Y)  EX  P  (- i  ETI  (J  X + m  Y))  A  X  d  Y 


(6.3) 


-40  —00 


Using  this  and  Equation  (3.8)  we  have  for  a  given  (£0,m0,f0)  that 

S  (KkS)  =  &***(*-  ci 

E*P  (-  i  IV  (1. X  -*■  K>Y))  Jx  (i  V 

E*p(-iEir(x(J«-l)+Y  (Tno-^AxA^cDim 

.M  . ak  *(M 

«o  « 

=  EX  P^-i  £1t  (x  (i-S) + Y  Jiw 

-  ao  -»  *" 

—  —  (6.4) 


As  expected,  S(£,m,f)  is  a  two-dimensional  convolution  of  the  true 
directional  spectrum  S(£,m^f)  with  G(£,m)  the  Fourier  transform  of 
g(X.Y).  We  see  then  that  S(£,m,f)  is  a  weighted  average  of  S(£,m,f) 
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and  that  G(£,m)  is  a  measure  of  the  directional  resolving  power  of  the 
assumed  probe  array.  By  the  nature  of  g(X,Y)  we  have  from  Equation 
(6.3) 


G  (!,>»)=  1  +  ZfL  £  c 05  ( E Tt  ( l  M  X j)) 


(6.5) 


If  we  assume  that  the  wave  equation  (Equation  (3.1))  holds,  we  find 
that  S(£,m,f0)  is  zero  when  l  +  m  ^  k0,  the  wave  number  for  fQ. 
From  this  we  have  for  a  given  (£0m0)  that  the  directional  resolving 
power,  DRP,  is 


M  CO 

DRp(l,*i,f.)  =  G  A  <A>n 

ter  J  =:  K0  Cos  0  ,  hi  =.  K.  SIN  0, which  implies 

and  we  get,  for  energy  coming  from  a  direction  6Q  at  frequency  fQ  as  a 
function  of  0  <_  0  _<  2tt  ,  that 

DRP  (e  1 00, f.)  =  G(k.(cos  0o-cose)3K.(siN0.rSiN0))  . 

From  Equation  (6.5)  we  get 

DRP  (e  |e.,(.)=  1 + 2*£  Xcos[ati(xiJK.  (cos  e. -cos  e>  K. 

i  =  J  jcxi-i 

,  ,  (sin  00-  -Sid  £))] 

=  1  +  E^^jcosklT ((i«-K.cos  (m.-  K.  Sin 


note  that 


*  =  l  j-A+l 

HAT  DRP(9,|  90fJ  =  N  (N  - 1)  +  I 


(6.6) 


Compare  Equation  (6.6)  and  (5.5).  Except  for  the  amplitude  term,  A^/4M, 
in  Equation  (5.5)  the  equations  give  identical  results.  The  choice  of 
bij  =  [N(N-l)  +  l]“b  =  1/M  is  again  found  convenient. 
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7.  DIRECTIONAL  ANALYSIS  FROM  THE  CROSS  SPECTRAL  MATRIX 


First,  we  consider  a  fundamental  approach.  We  have  for  a  pair  of 
detectors  I  and  J  located  at  (x^.y^)  and  (x^yj)  respectively,  the 
cross  spectral  matrix,  P^j (f)  =  C^j (f)  +  iC^ . (f)  and  more  impor¬ 
tantly  4>ii  (f )  the  phase  lead  of  I  over  J  given  by, 


<P.(f)  =  AB.tr AH  |g£lj 


The  actual  phase  lead  may  differ  from  this  value  since  the  true  phase 
lead  0  is  some  one  of  the  values 

%(i)  +  fi  21T 

where 


Consider  a  single  wave  of  frequency  fQ  with  corresponding  wave 
length  XQ  and  wave  number  K0  =  1/XQ,  and  find  the  direction  the  wave 
must  travel;  i.e.,  fit  a  single  wave  to  the  spectral  matrix  results  for 
the  detectors  I  and  J. 

Let  Dj_j  be  the  distance  between  I  and  J.  The  distance  between 
I  and  J  in  wave  lengths  is  K0D.jj .  In  radians  this  is  2nK0D^j .  From 
this  relation  we  get 

-ZirK.Dy*  <M  K.Dii 

or  that  only  values  of  h  such  that 

-  ett  K.D (PJf)  +  K  atr  ^  zn  k.d4j  (7.2) 

give  physically  acceptable  candidates  for  the  value  of  $.  If  D-y  <  XQ/2 
only  one  h  value  is  valid.  If  AQ/2  <  D^j  <  XQ  at  most  two  values  of  h 
are  valid,  etc. 

The  problem  is  how  to  find  the  true  direction,  0O,  of  a  single 
wave  given  the  above  possible  value  (s)  of  <j>.  Consider  a  given  value 
of  4>  in  terms  of  wave  length  units  and  we  obtain 
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Figure  4  illustrates  a  case  for  $  >_  0  (and  thus  L0  0) .  From  the 
figure  we  have,  where  the  true  direction  is  0O  and  true  phase  is  <p, 
that 


Oo 

MHCRC 

OR 


UI..+  3L 

T«j  £ 

Sin  o L 


+  oC 

=  _i_ 

ZT1  K, 

r 


l 


OC  =•  ARC  SIN 


ARCS  1 


Recall  that  sin(a)  =  sin(+TT-a)  .  Thus  for  a  given  <f>  _>  0  we  get,  since  a 
must  be  obtained  from  an  arcsin  relationship,  two  possible  values  of 
wave  direction,  the  true  value  6Q  and  its  image 


«'=[to+f]  + 

=  to-  i  -«•  =  Uu-[$  +ot] 


This  is  illustrated  in  Figure  5.  In  actual  practice  we  do  not  know  the 
true  direction,  thus  a  given  value  of  <j>  >_  0  gives  the  direction  as 

e.=  HI,*  [f +eC] 

where  0  <_  a  _<  it/2  is  obtained  from  the  principle  value  of  the  arcsin. 

If  (f>  <  0,  then  I  actually  lags  J  by  |  <j>  |  >0  and  Lq  =  Lq/D^j  <  0,  so 
that  arcsin  (Lq)  gives  -tt/2  <_  a  <  0. 
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The  implied  directions,  for  <j>  <  0,  are  directly  opposite  from  those 
for  <j>  >_  0  so  that  for  a  given  value  of  <j>  <_  0  we  get  directions  (a  <  0) 


e.  =  f  - 


or  A  —  U)L  .  +  3L  + 

°  Z 

6  =  >^*[1  + 

as  before. 

Thus,  where  the  principle  value  of  arcsin  is  assumed,  we  get  a  set  of 
possible  directions 


a  =  Yij  *  [f  +  SI 


where 


06,  =  ARCSIN  ' 

k  [_  air  K0Dij 

h  being  constrained  by 


(7.6) 


I  **  I 

An  example  of  this  type  analysis,  for  an  array  pair,  is  illustrated 
in  Figure  6.  Thus  several  estimates  of  0O  are  available  (at  least  two) . 

The  estimates  of  true  direction,  0O,  often  vary  from  one  array 
element  pair  to  another,  making  the  selection  of  a  true  0O  value  diffi¬ 
cult.  The  selection  is  also  hindered  because  half  of  the  estimates  of 
0O  are  of  the  image  type;  i.e.,  false  estimates. 

While  the  above  directional  method  leaves  something  to  be  desired, 
it  does  illustrate  the  basic  directional  information  produced  by  an 
array  of  detectors. 

A  better  method,  suggested  in  Munk  et  al  (April  1963),  of  using  the 
spectral  matrix  directional  information  to  fit  a  single  wave  at  each 
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FIGURE  6.  DIRECTIONAL  ESTIMATES  FOR  A  PAIR  OF  ARRAY  ELEMENTS 


frequency  is  given  below.  It  is  based  on  Equation  (4.16)  in  the 
form 

=  Jrfc.(fW2££rCii(0cos(zn(J,Kij+wY4 

r\  L  istirUiL 


-  Q.  ({)  sin  (a/tt  ( A  Xjj+  >1 '(;)')]] 


-D  +  1],  C.(4)  =  £  £  c,m 


where  M  =  [N(N 


and  N  =  number  of  array  elements . 


Recall  that  for  a  single  real  wave  of  frequency  fQ  >  0  and  known 
direction  60,  Ci;j(f0),  .  (fQ)  are  known  (see  Equation  (5.3)).  These 

values  give  J 

S/K,.e0.O  =  =  -£  • 

When  a  single,  well-directed  swell  is  expected,  it  is  reasonable  to 
assume  a  single  wave  for  a  given  frequency  exists,  and  to  select  0O  and 
A0  =  A^/4  such  that  the  least  square  error  between  the  theoretical 
cross  spectral  matrix  for  a  single  wave  (see  Equations  (5.1)  and  (5.3)) 
and  an  observed  cross  spectral  matrix  is  a  minimum.  Accordingly,  using 
the  expressions  of  Equation  (5.3)  and  an  observed  cross  spectral  matrix 
for  a  given  frequency,  fQ,  we  can  form  the  squared  error 

H  -  <C„  -  A„)2 


♦  z,££fCjj-A.cos(afl(i.  V 

+  a  £  £.  r<V  A.  sin  ( an  ( i.  X(i+  ».  Ya))l 

isijsi+iL 

\  N 

note:  A0=  ~  :  C0  =  -jjj-S  C.. 


(7.8) 


Expanding  and  collecting  terms  we  get 

H  =  Cp,+  2.22^1  jCi j  +Qij)  +•  +  lj  Ao 

-  iA.Jco  +  Cii  cos (iv  (i0  X.> 

"  §,VIN  V 

or  using  results  in  Equation  (7.7) 


(7.9) 
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A-l4“Ut  **  **K 


To  find  A  that  minimizes  H*  consider 

o  .  ■  • 


iii  =.  £[n(n-iV  i][Ao  -  =  o 


(7.11) 


This  requires  that  AQ  be  of  the  form  A0  =  S(L,m,fo)  and  a  resulting 
value  of  H  of  the  form 


k = c;  +  2. 


+q')-[n(n-i)+i][S(1,v,^  fjj* 


(7.12) 


2  2  2 

Since  CQ,  C^j  ,  and  are  all  nonnegative,  a  minimum  H  results  when 
S(S,,m,f0)  is  a  maximum.  A  choice  of  lQ  and  mD  that  maximizes  S(£,,m,f0) 
implies  a  0O  =  arctan  m0/$,0  which  is  optimum.  Remember  that  we  are 
assuming  holds,  along  with  the  wave  equation.  The  results 

then  for  each  fQ  is  a  two-sided  energy  spectrum  estimate  A  (fo,0o). 
Appendix  B  contains  a  listing  of  a  FORTRAN  II  program  for  finding 
Ao(fo,0o),  the  least  square  wave  fit,  from  a  set  of  spectral  matrices 
obtained  from  the  task  SWOC  data  collection  and  analysis  system 
described  in  Bennett  et  al  (June  1964) . 

Examples  of  least  square  single-wave  train  fit  analysis  from 
Bennett  (March  1968)  are  shown  in  Figure  7. 

A  more  complete  collection  of  the  directional  spectra  calculated 
from  data  collected  off  Panama  City,  Florida,  is  given  in  Appendix  A, 
and  in  Bennett  (November  1967),  and  Bennett  and  Austin  (September  1968), 
an  unpublished  Laboratory  Technical  Note  TN160. 

We  would  actually  like  a  continuous 'estimate  of  S*(Jl,m,f).  Con¬ 
sider  then  a  third  method.  From  Equation  (3v8)  we  have  for  a  pair  of 
detectors,  as  illustrated  in  Figure  8,  that 

a  a  ~ 

P*(X,Y,4)  =  [  fs(l.>..{)Exp[i2*(JX  +  YnY))  <U<U  . 
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DATE:  I  SEPT.  IMS 
TIME:  IS29-ISS2 


FIGURE  -^DIRECTIONAL  SPECTRUM  Ao(f,0) 


DATE:  4  JUNE  IMS 
TIME:  1411-1449 


FIGURE  7<.  DIRECTIONAL  SPECTRUM  A0(f ,  91 


DATE:  4  JUNE  IMS 
TtNE:  ISM- ISSt 


FIGURE  DIRECTIONAL  SPECTRUM  A0(  1, 9) 


DATE:  9  SEPT.  IMS 
TIME:  1129-1153 
STAGE  II 


FIGURE  70  DIRECTIONAL  SPECTRUM  A0(f,9) 


FIGURE  7.  LEAST  SQUARE  SINGLE  WAVE  FIT  DIRECTIONAL  SPECTRA  A  (f,6  ) 

o  o 
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Using  the  above  changes  of  variable,  we  get 

P*(x,Y,0  =  J  [s(K  Q  EXpfian  (KDcosecosiJj+KosiNesm 
m1t  KAkde 

If  «o 

P  X  ,Y, ^  =  j  (K,e,f  )  E  X  p[i  air  K  D  c  OS  (e  -  (p  | K  A  K  d  0 


We  have  from  Equation  (3.2)  that 

£(K,e,f.)  =  S(k, 

If  we  think  in  terms  of  fQ  >  0 

^  =  a  S  (K,0,fo)  (7.i4) 

We  are  assuming  that  the  wave  number  relation  of  Equation  (3.1)  holds. 
Thus  Figure  3  is  applicable,  and  we  can  write  the  one-sided  spectral 
density  as 

sVk.6,0  =  £  <£  (k-k„) 

i(K-K.)={'  K=K- 

Co  **  K« 

This  allows  us  to  write,  — .  ^  ^  ^  , 


teJHtftC. 


•Tf  CO 


r>;u) =]{*(«.  f.)«f(K  -K,)Enp[i.2-rr  K  Dcos(e  -  iji^KolK  Afl 


-It  0 
* 


P*(xxU=  [a(e,f.)  txp[i  t  it  K„Dcosfe-ii))jk„le 


-Tf 


(7.15) 


We  have  reduced  the  problem  to  finding  [a(0,fo  •  KQ] . 

From  Equation  (7.15)  we  see  that 

?*(o,o,S\-  P(f)=  ja(©,t.)K.Jie 

-If 

where  P(f0)  is  the  power  spectral  density  of  frequency  fQ.  For  better 
comparison  of  cases  where  P(f1)  =  P(fQ)>  jf-J  4  [fQ|,  it  is  convenient 
to  express  a(0,fo)  in  a  normalized  form 

A(0,fo)  =  a(0,fo)Ko 
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where  we  get 


IT 

p(£>  '  (AM)Ae  (7.i6) 

-nr 

Thus  if  the  energy  distribution  as  a  function  of  direction  is  the  same 
for  P(fl)  =  P(fQ)  then  we  also  get 

ACe.fi)  =  A(0,fo). 

Consider  now  (assuming  A(0,f)  can  be  so  expressed)  a  Fourier  series 
expansion  of  A(0,f)  for  fixed  f.  Clearly  it  is  periodic  in  0  with 
period  2tt.  Thus  for  any  given  f  =  fQ  we  can  write  A(0,f)  as 


oo 

t\(e)  =  ~  ne  t  kslN 


(7.17) 


Substituting  this  expansion  into  Equation  (7.15)  we  get 

r  ir 

?*(X,Y,0  =  EKP  i  2-tr  KD  C0s(tf-^)  Ad  +  [cos  Y\9 

Cm  J  H*l  J 

-HT  «  -| 

Ekp  i  2.  It  Kfc  cos  (d-ty)  cU  +  by,  JSU4  me  E*P  lair  KD  cos/d-Yji® J 
or  P*(X,X*)  =^fjcos(aitKD  cos(e-i|>))oU  COS  >10  CoS  (2tC  KD 

cos(e-\jf))Ae  cos(ETcKDcos(e-i|>)yie)] 

f  ■lr 

+  AlT  r‘>l  KD  C0Std-t))cL®  + 

-tr 

¥yl&Ato$yies\*i(WK*Co$ie-iifa  + 

•«=i  '  Mr 

b„jsi»i  ne  Sin (a|Ko  cos(#-i|r^Ae^J 


(7.18) 
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Thus  we  can  express  P*(X,Y,f )  as  complex  infinite  series  with  unknown 
coefficients  3q,  a*,  a2»...  and  b^,  b2>  ...  and  constants  defined 
by  the  integrals  (let  Z  =  2ttKD;  n  =  0,  1,  2,  ...)  of  the  form 


^  cos  6  c os  (i  Cos  (0  -^))  Ae 
^SlM7|0  CflSUCOS^-Tlr^Ae 
Jcos7i0SiM(fLCos(0-7jr5)Ae 


aid 

i-tr 


[sin  no  SlM(aCOS(0-^A0 

'•K 


(7.19) 


(7.20) 


(7.21) 


(7.22) 


Consider  Equation  (7.19)  where  $  =  (8-if0  and  d<f>  =  d6,  ip  being  a  constant. 
We  then  have  on  changing  variables 


(7.23) 


Since  the  integrands  in  Equation  (7.23)  are  both  of  period  2n  and  the 
interval  [-tt-i/;,  is  of  length  2n  we  can  write  the  equivalent  of 

Equation  (7.23)  as 


COS  VhMCOS  Y)(f  cos(i  C0$^)d$ 

-71  If 

—  Sim  Isin  yi4cos(zcos^)cl(p  (7.24) 
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From  Ryzhik  and  Gradshteyn  (1965  -  page  402)  and  noting  that  the  second 
integrand  is  odd  we  get  Equation  (7.23)  equivalent  to 


COS  Y[ 


(7.25) 


where  Jn(Z)  is  the  Bessel  function  of  the  first  kind.  Employing  a  simi¬ 
lar  procedure  for  Equations  (7.20),  (7.21),  and  (7.22)  we  get 

P*(x,y,0  =  J.(ankD)j  + 

£ [a* cos  rity  sir  kd)  +  b„si*  ny  zit 

+  iZ[a’<C05,>'l'  J„(2hkd) 

+  bn  Sin  n  iy  Zir  Sin^OI^  Jn(2H  K  D)1 


(7.26) 


Now 


and 


cos/20i\  -  J  0  * 

UJ  -  U-0%  n 


oddL 
Cv/CV> 


$\M 


Thus  we  get 


/r?Tt\  _  f  0  n  even 
\  2  /  |  (-()¥  n  oAA 

P*(X.Y.O  =  §*[»  J.UUKB')] 


+  £  [in  J„(tH  kd)  (- 1)*  (a„C0S  V|l|l  +  b„  sm  * ljf)l 

*1*  1,4,6,— 

+j£[z-nJn(iTtK0)(-i)'‘jJ  (anC»5  ?iy+^siNioiv)j 


ii*  *.M- 


(7.27) 
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From  a  spectral  matrix  of  the  form  in  Equation  (4.13),  M  =  N(N-l)  +1 
different  equations  can  be  set  up  using  Equation  (7.27).  This  allows 
us  to  get  a  system  of  equations  for  any  m  of  the  unknown  coefficients 
aQ,  a^,  a2>  •••  ;  b^,  b2»  b3,  ...  while  assuming  the  rest  of  the 
coefficients  are  negligible.  We  can  then  solve  for  the  m  desired 
coefficient  values.  This  has  not  worked  well  in  practice  for  two 
reasons.  The  inverse  of  the  matrix  of  constants  obtained  is  sparse 
and  often  ill-conditioned.  Further,  if  the  wave  energy  is  from  a  nar¬ 
row  beam  width  (30  degrees  or  less),  the  first  100  harmonics  in  the 
Fourier  series  expansion  can  be  significant.  There  is  perhaps  a 
more  efficient  orthogonal  set  of  functions  than  the  sines  and  cosines 
of  the  standard  Fourier  expansion.  Search  for  such  an  orthogonal  set 
should  prove  fruitful.  One  might  start  with  Walsh  or  Haar  functions. 
See  Hammond  and  Johnson  (February  1960) . 


8.  SUMMARY 

It  is  believed  that  the  least  square  method  of  using  the  informa¬ 
tion  in  a  spectral  matrix  is  the  best  method  presently  available. 
Examples  of  such  analysis  can  be  found  in  several  of  t)ie  papers  in  the 
bibliography.  A  collection  of  ocean-wave  induced,  bottom  pressure 
directional  spectra  from  these  papers  is  given  in  Appendix  A. 

An  iterative  extension  of  the  least  square  method  can  be  found  in 
an  excellent  paper  by  Munk  et  al  (April  1963) .  Some  details  of  this 
method  are  given  in  Appendix  C  along  with  an  example  result  and  a 
FORTRAN  program  for  the  method. 

There  is  merit  to  using  the  coherency, 


to  form  the  weights  b. .  in  Equation  (4.16).  One  idea  being  explored  is 


'  tj-i  fJ 
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APPENDIX  A 


A  COLLECTION  OF  DIRECTIONAL  OCEAN  WAVE 
BOTTOM  PRESSURE  POWER  SPECTRA 


This  appendix  is  a  collection  of  the  results  of  a  least  square, 
directional,  single-wave  train  analysis. of  the  cross  power  spectral 
matrix  resulting  from  the  analysis  of  ocean  bottom  pressure  data.  The 
data  were  collected  at  Stages  I  and  II  offshore  from  Panama  City, 
Florida,  during  1965.  The  data  collection  system  and  the  estimation  of 
the  cross  power  spectral  matrix  associated  with  a  set  of  data  are 
described  in  Bennett,  et  al  (June  1964).  Augmented  pentagonal  arrays, 
containing  six  pressure  transducers  each,  were  located  seaward  of  each 
of  the  stages.  Stage  I  is  11  miles  offshore  in  approximately  103  feet 
of  water,  and  Stage  II  is  2  miles  offshore  in  approximately  63  feet  of 
water. 

Certain  parameter  values  are  pertinent  to  the  directional  analyses 
presented:  the  number  of  data  points  in  each  pressure  data  set  is 

N  =  1800;  the  sampling  rate  is  once  per  second.  At  =  1  second.  On  the 
cylindrical  polar  plots  frequency  is  the  radial  variable  and  compass 
bearing  the  angular  variable.  The  vertical  axis  is  log^g  of  power 
spectral  density  in  inches  -seconds  of  water  pressure.  The  frequency 
axis  range  is  0  to  0.3  Hz  in  0.05  Hz  increments.  This  is  illustrated 
in  Figure  7  of  the  report.  In  each  plot  title,  the  date,  time,  and 
location  (stage)  is  indicated.  The  value  WD  is  wind  direction  in  com¬ 
pass  degrees,  and  WS  is  wind  speed  in  knots.  Appendix  B  gives  a  listing 
of  the  FORTRAN  II  computer  program  used  to  produce  the  plots. 
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DATEO6/04/65  TIMEI4I8-I449  STAGEIWo282 
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DATE  06/05/65  TJITO33-OIOI  STBGtl  WIOO 


OATE06/05/65  TJ»CD  1 02-0130  9TMZI  MO  1 00 
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APPENDIX  B 


A  FORTRAN  II  PROGRAM  FOR  SINGLE-WAVE  TRAIN  ANALYSIS 


The  FORTRAN  II  listing  of  an  IBM  704  program  for  the  least  square 
single-wave  train  analysis  of  a  spectral  matrix  obtained  from  an  array 
of  ocean  wave  bottom  pressure  transducers  is  included.  The  listing  is 
from  the  FORTRAN  to  ALGOL  translator  of  the  Burroughs  B-5500  and  is 
syntax  free  at  the  FORTRAN  II  level.  The  mathematics  of  the  single¬ 
wave  train  analysis  is  described  in  the  body  of  this  report.  The  data 
collection  system  and  calculation  of  the  required  cross  spectral  matrix 
is  described  in  Bennett  et  al  (June  1964).  The  plotter  subroutines 
GPHPVW  and  PFB3D  are  also  included.  Bennett  (March  1968)  describes  the 
plotting  technique  in  PFB3D  which  was  used  to  produce  the  plots  in 
Appendix  A. 


FORTRAN  TO  ALGOL  TRANSLATOR 
PHASE  1  FORTRAN  STATEMENTS 

C  LISTING  OF  AN  IBM  7 OA  FORTRAN  PROGRAM  FOR  DIRECTIONAL  WAVE  ANALYSIS. 

C  13350-2  C  BENNETT  JULY  1967  SWOC  78*301-8210 

c  Request-0268 

C  LEAST  SQUARL  SINGLE  WAVE  TRAIN  FIT  AFTER  MUNK 

DIMENSION  IOC  11 >*FQ( 100)»WNt 100 ) *  PER  I 0D( 1 00 ) > W A VLGH ( 1 00 ) > THM A  X ( 100 
X)»DEpATN( 100)/ AVEP(100)#P(100»6»6)#EMAX(100)»TEM(146) 

DIMENSION  D(6,6>.PS!(6,6),E(100,72),THETA(72)/DE(73>#THZRU(1A6), 
XH(100)»Hl(l00)#H?(10O)  / GOOD{ 6 ) 

DIMENSION  data  (1000) 

E2TEN=0 , 43429448 
TW0PI*6. 281853 
RT0*S7. 29578 
DTR=0. 0174532925 
FIV*0, 087266465 
REWIND  3 
REWIND  9 

CALL  PLOTS(DATA(tOOO>MOOO) 

CALL  PL0T(0.0#-30,0>-3) 

CALL  PL0T(2.S#2.5/-3) 

DO  5  1*1  >72 
Zl»(l-1)*5 
5  THET  A( I )*Zl *DTR 
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00  10  1*1*6 
DO  10  J*1 *  6 
0 ( I »  J  )=0  .  o 
10  PSI(I* J3=0.0 
0(1*23=100,0 
0(1*33=100.0 
0(1*43=100.0 
0(  1*53*100,0 
0( t  »  6  3  =  1  00 ,0 
0(2*33=117,558 
0(2*43*190.212 
DC  2  .“53  =  1  90.21  2 
0(2, 63=11?. 559 
0(3*43=117.559 
0(3*53=190.212 
0( 3*63=190.212 
0(4*53=117.558 
0(4.63=190.212 
0(5*83=117,559 
C  PSI  IS  TRIG  ANGLE 

C  0  IS  DISTANCE 

pSI( 1*23=0.0 
PSI(1*33=72.0*DTR 
pSI ( 1 *43=144. 0*0TR 
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f,SI(l>5)  =  2l6.0*DTrt 
PSICl»6)=2e0.O*DTR 
PSI(2»3)-126.0*0TR; 

PSIC2>4)=162.0*DTR 
PSIC2»5)=196.0*0TR 
PSI(?>6  )=234.0*0TR 
PSI( 3»4)=196.0*DTR 
PSI(3»5)=234.O*0TR 
PSI( 3>6)=270.*0TR 
PSIC A*5)=270.0*DTR 
PSI(4»6)=306.0*DTR 
PSlC5»6)  =  3fl2.0*OTR 

C  GOODCl)=  1  FOR  USARLE  channel  data  AND  o  for  BAD  CHANNEL 
15  READ  2,  ISK1P  » ( GOODC I )* I31 >  6 ) 

2  format ( i 2  ,6Fi.o  ) 

IFCISKIP)100»30»20 

20  DO  25  I*1MSKIP 

READ  TAPE  3»ID»FQ(1 )»WN(1)»M>K 
IF (K)2l»22»21 

21  PAUSE  20202 

C  TAPE  Out  of  phase  WITH  data  read  DESIRED 
GO  TO  15 

22  MP=M+1 

DO  23  L  =  2  »  MP 
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23  READ  TAPE  3 
25  CONTINUE 

fill  TO  15 

C  F0(l)=0.0  CAN  NOT  BE  ME  A N I NflL YE UL L Y  PROCESSED 

30  RE  AO  TAPE  3»I0,FQ(U*WN(1),M.K 
I  F  (  K  >  2  1  *  3. 1  >  2  1 

3 1  MP  =  M+1 

no  32  L  =  2»  MP 

32  READ  Tape  3,ID,FGCl).RNCL),M,K,DEITaT,i;>EPTH,PERI0D(L),WAVlGH(L), 
XDEPATN(L)»AVEP(L)»((P(L>I>J)»J=W6)>I  =  i>6) 

C  P(LM'J)=*0,0  IF  CHANNEL  1  OR  J  IS  NO  GOOD 

c  wn(L)  is  wave  mjmbeb=2pi/w ave  length  in  feet 

C  PI  =  3. 141592/.  .  . 

C  VALUE  K  NO  LONGER  neeoed  , 

GOODCH=0.0 

I 

DO  AO  I»l>6 

G000(  I  )*=GOOD(  I  )*P(2*  I  >  I  ) 

IF(GOOD( I)  ) A  1 . AO. 41 
A  1  GOOOCH  =  GOOOCHM  ,0 
GOOO ( I )  =  1 . 0 
AO  CONTINUE 

terms=i.omgoodch*cgoodch-i.o')  ) 

00  70  L=2»MP 
SUM=0.0 
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DU  50  I  =  1 » 6 
I F  (  GOOD  (  I ) )  b  1  .50.51. 

51  5UM  =  SUM  +  PCL.  I  .  I) 

50  CONTINUE 

AVEP(L)=SUM/GOODCH 
DU  60  K  =  1 »  72 
SUMsO.O 
DO  90  1=  1*5 
IFCGOUD(I)  )90»  90»  91 

91  IP  *  IM 

DO  95  J«  I P  *  6 

IF(GOUd(J)  )95,95»q? 

92  SUM  *  SUM  +  P(L»I#  J)*C0sFCWN(D*DC  I#  J)*C0SF(THETA(K)-PSI(  I»  J))  ) 

X  -Pel* J#I)*SINF(HN(L)*D(I#J)*COSFCTHETA(K)-PSi(I# J>>) 

95  CONTINUE 
90  CONTINUE 

60  E(L#K)=( AVeP(L)+2.0*SUM)/TERMS 

i f ( sense  switch  U62.63 

62  PRINT6,FQ(L) » AVEP(L1 • ( I D ( I) # I  *  1# 1 U , (E (L#K ) » K»  l  , 72 ) 

6  FORM aT(1Xm PEI  1.AMPE12.4M1A6/(10(1XMPE1 1.4))) 

63  G=0.1/0TR 

C  DELTA  THETAe5DEG  UR  H*1/2*DEL  THETA 

DE(1)°G*(E(L»2)*E(L#T2)) 

DEC72)=G*(E(L.l )-E(L»71  ) ) 
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DO  6 1  ITHETA=2*71 

61  0£(ITHETA)eG*(E(1*ITHETA+1>-E(L»ITHETA-1>) 
DEC  73  )  =  DE(  1 > 

K  =  1 

DO  65  1=1*72 

IFCDE(I)*OE (1+1)566*67*65 

66  71=1-1 

C  FIV  IS  5  DEG  IN  radians 

THZR0CK)  =  FIV*CZI  +  CQFC  n/CDEC I) -DEC I  +  l ) ) ) ) 

K  =  K  + 1 

•  GO  TO  65 

67  IFCOFCI >>60,69,60 

69  THZROCK)=Fl V*FLOATFC 1-1  ) 

K  =  K  +  1 

60  IFCDECI+l >>65*64,65 
6a  THZROCK)  =  FlV^FLOATF(  I) 

K  =  K  + 1 

65  CONTINUE 
NZEROS*K-l 
DO  75  K=l»NZEROS 
C  TERMS  SAME  AS  ABOVE 

SUM»0 • 0 
DO  76  1=1*5 


IFCGOODCI)  >76*76*77 


77  I P  =  I ♦  1 

no  79  J=  I P  *  6 
IFCGGOOCJ)  )79*70*78 

7 8  SUM  =  SUM*P{L,I,J)*cnSF(WN(L)*0(I»J)*COSPCTHZRO(K)-PSICI#J))) 

X  -P(L. J» I )*SINF(WN(L)*D( I> J)*C0SFCTHZR0(K)-PSI(I» J))) 

79  CONTINUE 
76  CONTINUE 

75  TeM(<)=cAVeP(L)+5.o*SUM)/TERMS 
EMAXC  U  =  TEM(  1  ) 

THMA  <Cl)=RTU*C  TWnPl-THZROCt  )  ) 

DO  7A  K=2>NZERGS 
IFCEMAX(L)-TE^(K))73»74»7A 
73  EMAXCL)=TEMCKT 

thmaxcl)=Rtd*ctwopi-thzrO(K)  ) 
c  thmax  is  bearing  from  magnitic  north 

7U  CONTINUE 
SUM=0.0 
00  80  1=1*5 
IFCGOUDCI)  )80*80»81 

81  IP=I+1 

DO  83  J  =  I P  »  6 
IFCGOOD(J)  )83»83»82 

82  5UM=SUM+P(L* I* J)*P(L* I* J)+P(L* J* I )*P(L* J* I) 

83  CONTINUE 


B-8 


80  CONTINUE 

SS0»AVEP!L)*AVEPCLU2.0*SUM 
H(LjaSSO-TERHS*EMAX(L J*EMAX(L) 

atilda»av£pcu/twopi 

HTIL0A«S$0-ATILDA*ATILDA*TERMS 

Hl(L)=1.0-(h(L)/NTlL0A) 

H2CL)*CEMAX(L)-ATILDA)/(AVEPCL)-ATILDA) 

if( sense  snitch  i)99»ro 
99  PRINT?’#  TERMS.  EMAX(L)#ATlLOA#HTlLi)A 
r  formaT(ix#f5.i#ip3e:ii.a  j 
70  CONTINUE 

T  OF  # 1 1  A 6 / 

» 1 IHATTENUaTION. 4X» 
2HH2) 

EP(L>  #EMAX(L). 

•  A# 

0PF6.3.6X.0PF6.3  ) 

00  1 l 0  U=2#MP 

AVEPCU)=E2TEN*L00F{A8SF(AvEPCU)  ) 

EMAXCLjaE2TEN*L0GF(ABSFCEMAXCL))) 

UO  WAVLr,H(L)=EMAX(L)-E2TEN*L0GFC  ABSF(DEPATNCL))) 


PRINT  3# i ID!  I )»  I  =  l#l  1  )  »!G000! I)»I=1»6).M 

3  FORMAT! 1H1,39HLEAST  SGUARE  SINGLE  NAVE  TRAIN  FI 
X1X#6F2,0*2X#2HM=I3// 

X3X#9HFREQUenCY»5X#6HPERI0D»3X.1 IHWaVE  LENGTH* 1 X 
X5HAVE  P#8X,1HA,10X,5H8RNG  #8X» IHH# 1 2X #2HH 1 , 1  OX# 
PRINT  4#(F0cL)#PERi00CL)»WAVLGH(L)#0EPATN(L)» av 
XTHMAXCL)»H(L  J»H1(U#H2(L)#L*2»MP) 

4  FORMAT! lX#0PFlt.T»3X#0PF9. 5# 1X»0PF11. A#1X#1PE11 
X  lX.lPEll.A*  !X#lPEll  .4#3X,0PF7.2#3X»lPEn.«»«X, 
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AVEP( 1)=AVEP(2) 

EMAXC 1)=EMAXC2) 

WAVLGrt(l)=KAVLGH(2) 

CALL  SyMB0L(2. 0,0. 75*0. 1,14HWAVE  TRAIN. FN=  ,0.0,1a) 

CALL  NUMBER(3.5,0.75»0.1»FQCMP),0.0,2) 

CALL  GPHPVW(MP,in,AVEP*EMAX,WAVLGH) 

c  wavlgh  is  surface  wave  staff  spectrum  EST  FROM  bottom  pressure. 

C  AvE=P  EMAx=V  WAVLGH=W  ON  THE  PLOTS. 

ZM  =  M 

NP=DELTAT*ZM*0.5+1 .0 

DELX=0.0 

DO  115  1=1,6 

CALL  NUMBER ( DELX,0. 1 ,0. 1 ,GOODC I  )  ,0.0,-l  ) 

115  DELX  =  l)ELX  +  0. 1 

CALL  PFB3DCNP,I0.EmaX,FQ.THMAX> 

GO  TO  15 
100  REWIND  3 
CALL  EXIT 
PAUSE  70707 
GO  TO  15 

C  .  EN0(0»1, 1,0,1  ) 

END 


B- 10 


FORTRAN  TO  ALGOL  TRANSLATOR 
PHASE  1  FORTRAN  STATEMENTS 

SUBROUTINE  GPHPVWfL, I D* ZLOGP , ZLOG V , ZL OGW 1 
C  SPECIAL  FORM  OF  GPHPVH  FOR  13350-2  AUGUST  1967 

n I  MENS  I  ON  ZLOGP(?)*ZLOGV<21,ZLOGW( 2 1*10(21 
CALL  PLOT(0. 0*0. 0*31 

call  plot  (0.0*10.5,21 

call  plot  (8.0*10.5*21 

call  plot  ce. 0*0.0,21 

call  plot  (0.0*0.0*21 

call  plot  (2. 0*1. 5, -3) 

call  symbol  co.o.-t  .0,-.  1#  io(ii*o.o*66i 

call  AXIS  (0.0,0.0,20HNORMALI^EO  frequence .-20,5.0.0.0*0.0.0.21 
7MAX  =  ZL0GP(  1 1 
DO  to  I =2*  L 

compak=zlogpc II 

10  ZMAX3MAX1F(ZMAX»C0MPAR1 
MAX=7MAX+3.0 
BE=MAX-8 

CALL  AXIS(0.0*0.0*17HLOG  POWER  DENS  I T Y *  1 7* 8 . 0* 90 • 0* BE , 1 . 0 1 
DX=5.0/FL0ATF(L-1 ) 

Y=ZL0GP(1 1-BE 
X  =  0.0 

call  plot c x * y * 3 i 
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00  20  I=2»L 
X*X+0X 

Y=ZLOGPC  I  )-6E 

20  CALL  PLOT  CX.Y.2) 

CALL  SYMB0L(X»Y»0.1»1HP*0.0»1) 
X  =  5.0 

Y=ZLOGVCL)-BE 

CALL  SYMBOLU.Y.O.t.lHA, 0.0.1  ) 
CALL  PL0TCX.Y.3) 

M-L“l 

DO  21  1=1* M 
X=X-OX 
1 1 =L"I 

Y=ZL0GV(II)-8E 

21  CALL  PLOT ( X#  Y  »  ? ) 

WMAX=HAX 

X  =  0.0 

Y  =  ZLOGW( 1  >-BE 
CALL  PLOT  C  X  »  Y . 3  > 

00  2?  I=2« L 

x=x+ox 

IFCZLOGWC I 5-HMAX 524.23.23 
23  Y=8 .0 
GO  TO  22 
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24  Y=ZLOGW(I)-BE 
22  CALL  PLOT ( X>  Y » ? ) 

CALL  SYMBOL(X»Y»0.1#1HS*0,0*1) 
IF  (SENSE  LIGHT  1 )40>30 
30  CALL  PLOT  ( -2 . 0 > 9 . 0, -3  ) 

SENSE  LIGHT  1 
GO  TO  50 

40  CALL  PLOT  ( 6 . 0* -t 2. 0»-3  ) 

50  RETURN 

C  END(0M*1»0,0) 

END 
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FORTRAN  TO  ALGOL  TRANSLATOR 
PRASE  1  FORTRAN  STATEMENTS 

subroutine  pfr30cn, io,p,r,B) 

C  P=FUNCTlON  OF  R  a  POWER  DENSITY 

C  B=  function  OF  R  s  COMPASS  bearing 
C  R=  FREQUENCY  IN  RZ 

DIMENSION  IDC2)*pC?)*R(2>*BC2),CC360)*SC360) 

T  =  1 0 . 0 

0=0.70710676 
0TH=0. 017453293 

C  DTH  IS  ONE  DEGREEE  IN  RADIANS  IE  RADIANS  PER  1  DEGREE 

A  =  0.0 

00  10  1=1*360 

As A+OTH 

CC I  )=COSF(A) 

10  Scn  =  SlNFCA> 

CALL  PLOTCO.O, 0,0*3) 

CALL  PLOTCO.Q* 10.5,?) 

CALL  PLOTCB.O. 10.5,?) 

CALL  PLOTCe. 0*0.0,?) 

CALL  PLOTCO.O, 0,0*2) 

CALL  SYM80LC1. 0*0. 5, -0.1*10(1), 0.0*66) 

CALL  PL0T(A,0,3.5,-3) 

CALL  PL0TC3. 0.0.0,?) 
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CALL  SYMB0L(3.25>0.0*0.1/1HN  zO.O/1) 

CALL  AXISCO.O/ 1 .O/OH  /  0 *  5 . 0* 90 . 0/ " 1 . 0/ 1 .0  ) 

CALL  PLOTCO.O/1 .0/3) 

CALL  PL0T(0.0>0.0/2) 

CALL  PLOT (-3. 0/0. 0/?) 

CALL  SYMBOlC-3. 26/0. 0/0. 1/1HS  >0.0/1) 

X=-3. 0*0 
Y  =  X 

CALL  SyMB0L(X-0.1/Y-O.1/O.1>1HE  >0.0>1) 

CALL  PL0T(X/Y>3) 

CALL  PL0Tt-X/-Y/9) 

CALL  SYMB0LCO.1-X/O.1-Y/0.1/1HW  >0.0/1) 

CALL  PL0T(0.0>0.0/-3) 

DO  20  1=1/6 
A  =  0.5*FL0ATF  C  I  ) 

CALL  PLOTCA/O.Oz 3) 

DO  30  J= 1/360 
Y=A*S( J)*D 
X=A*CCJ)  +  Y 
30  CALL  PL0T(X/Y>2) 

20  CONTINUE 

CALL  PL0TC0.0.0.O/-3) 

0=10.0*0 

C  1 0*D  NEEDED  TO  SCALE  0.0"0.3  TO  0-3  INCHES  ON  PLOTS 
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00  40  I  =2*  N 
RAD=DTh*B (  I  ) 

Y  =  R< I )*D*Sl NF{ -RAO  ) 

XeR(I )*T*COSf(RAO)  ♦  Y 
IUD=3 

CALL  PLOT(X,Y, IUD) 

AYaPCD  ♦  Y  +  1.06 

C  2  IS  not  added  sn  Tqp  SYMBOL  ^ 1 LL  BE  THE  REFERENCE  IN  .00  SYMBOL 
call  PL0T(X,ay,2) 

CALL  SYMB0L(X,AY,0.08»2>0.0.-2) 

CALL  PlOT(X,AY,3) 

40  CALL  PlOT(X,Y,?) 

I F C SENSE  LIGHT  1)41,42 
42  CALL  PLOTC-4. 0,7. 0,-3) 

SENSE  LIGHT  1 
GO  TO  50 

41  CALL  PLOTU.  0,-14. 0,-3) 

50  RETURN 

C  ENDC0M»1'0M  ) 

END 
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APPENDIX  C 


A  FORTRAN  PROGRAM  FOR  ITERATIVE  WAVE  TRAIN  ANALYSIS 


The  FORTRAN  listing  of  a  Burroughs  B-5500  program  for  the  itera¬ 
tive  least  square  multiple  wave  train  analysis  of  a  spectral  matrix 
is  presented.  The  mathematics  is  an  iterative  utilization  of  the 
single  wave  train  analysis  described  in  the  body  of  this  report. 

Following  the  single  wave  train  analysis  of  the  measured  spectral 
matrix  the  resulting  values  of  single  wave  train  power,  A,  and  wave 
bearing,  6,  are  used  to  find  the  spectral  matrix  that  would  occur  for 
such  a  wave.  Details  for  this  are  given  in  Section  5  of  this  report. 

From  the  above,  a  residual  measured  spectral  matrix  is  formed  by  sub¬ 
tracting  a  fractional  portion  of  the  single  wave  spectral  matrix  from 
the  previously  used  spectral  matrix.  The  residual  spectral  matrix  is 
then  single  wave  train  analyzed.  The  above  procedure  is  continued 
iteratively  until  a  specified  number  of  iterations  have  been  com¬ 
pleted  or  the'  residual  spectral  matrix  total  power  gets  smaller  than 
a  specified  value.  Table  Cl  shows  numerical  results  for  several  fre¬ 
quencies.  Figure  Cl  is  a  plot  of  the  bearing,  0,  and  the  ratio  of  A 
to  the  total  power  available  in  the  original  measure  spectral  matrix 
for  the  frequency  band  0.00833  to  0.24187  Hz.  The  iteration  parameters 
were  set  for  a  maximum  of  five  iterations,  a  residual  power  ratio  of 
0.1,  and  a  fractional  portion  value  of  0.1.  Both  results  are  for  data 
collected  for  task  SW0C  at  Stage  II  between  1220  and  1249  hours  on 
4  June  1965.  The  wind  speed  was  12  knots  from  a  bearing  of  280  degrees. 
In  Table  Cl  A,  H,  and  bearing  are  as  previously  defined.  AVE  P  is  the 
average  C^^(f)  for  the  residual  spectral  matrix  and  P  is  the  AVE  P  for 
the  first  iteration;  i.e.,  the  original  measured  average  power.  The  HI 
and  H2  values  are  measures  of  the  isotropicity  of  the  energy  represented 
by  the  residual  spectral  matrix.  HI  compares  the  least  square  value  of  H 
with  the  value  H1  that  would  have  been  obtained  if  the  wave  energy  were 
isotropic: 


HI  =  1  -  (H/H’)  . 

H2  compares  the  power  A  of  a  single  wave  fit  to  the  total  power,  AVE  P, 
and  the  value  A'  that  would  have  been  obtained  for  isotropic  energy: 

H2  =  (A-A' ) / (AVE  P-A') . 

Both  HI  and  H2  are  between  0  and  1,  the  lower  limit  is  for  the  isotropic 
case  and  the  upper  for  a  plane  wave  from  a  single  direction. 
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